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Abstract. - We present exact results for the spectrum of the fractional Laplacian in a bounded 
domain and apply them to First Passage Time (FPT) Statistics of Levy flights. We specifically show 
that the average is insufficient to describe the distribution of FPT, although it is the only quantity 
available in the existing literature. In particular, we show that the FPT distribution is not peaked 
around the average, and that knowledge of the whole distribution is necessary to describe this 
phenomenon. For this purpose, we provide an efficient method to calculate higher order cumulants 
and the whole distribution. 



Anomalous diffusion is a widely investigated phenomenon with an increasing number of 
applications in natural sciences [1-6]. Stochastic Levy processes serve as a paradigm for 
many unusual transport phenomena in which collective dynamics, extended heterogeneities 
and other sources of jumps with a long-range distribution play an important role and lead 
to anomalous diffusion. Levy flights are governed by rare yet extremely large jumps of 
diffusing particles. In the continuous limit, which interests us here, the Levy flight process 
is described by a fractional diffusion equation for the Levy flyer concentration field C{x,i), 
given by 

dtC(x,t) = -(-A) N C(x,t), 0<N<1, (1) 

where (— A) N is the Riesz-Feller derivative of fractional order 2JV, namely (— A) N = ^ x] -in 
[7]. Eq. (1) describes a diffusive process that is faster than normal diffusion (super-diffusive) 
[1], where the index JV characterizes the degree of fractality of the environment. For JV > 1 
there is no probabilistic interpretation of the equation as a diffusion equation, although for 
various values of JV Eq. (1) can have a different physical interpretation. For example, N = 2 
corresponds to the overdamped vibrations of a flexible rod, or to the out-of-equilibrium 
fluctuations of a slowly growing film in molecular beam epitaxy [8]. 

Eq. (1) has to be supplemented with appropriate Boundary Conditions (BC) encoding 
the properties of the boundaries, such as absorbing boundary conditions 

(-A)"C(±l,t) =0, 0<^<JV/2. (2) 

Note that this condition should hold for all real values of \i in the interval < [i < N/2. 
Although we will be interested here in the case < JV < 1, it turns out useful to consider 
this equation more generally for any JV, and to specialize to the desired range when drawing 
physical conclusions. 
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In the context of Levy nights, a classical quantity of interest in such systems is the Mean 
First Passage Time (MFPT) defined as the average time needed for a stochastically moving 
particle to reach one of the two absorbing boundaries x = ±1, when it was initially located 
at some point x m the interval. It can be shown [9] that the First Passage Time (FPT) 
distribution for the one dimensional bounded domain with absorbing BC is obtained as 

d f 

p(t\x ) = - — J C(x,t\x )dx, (3) 

where the notation includes x giving explicit reference to the initial condition C (x, 0\x ) — 
5 (x — x n ). In particular, moments of the distribution p (t\xo) are given by 

/•OO 

(t m )(x ) = / t m p(t\x )dt. (4) 
Jo 

The MFPT in the presence of two absorbing boundaries can be exactly calculated using 
Sonin inversion formula [10], resulting in 

< ( ><*°> = rwri)' (5) 

where T(x) is Eulcr's Gamma function. This problem has been further studied both ana- 
lytically and numerically; generalized to other boundary conditions [11,12], to semi-infinite 
domains [13], and even to complex media [14]. 

A general question regarding the average is whether it is also representative. A simple 
test is to compare the average to the standard deviation of the distribution, or more generally 
to higher order cumulants. However, so far no attention has been given to the variance and 
other higher order moments of the distribution p (t\x ) of FPT for Levy flights. Indeed, if the 
distribution is peaked around the average, then knowing the average gives a fare information 
about the process. This is exactly the case for the Gaussian distribution where the average 
and the variance summarizes completely the information about the distribution. However 
if the distribution is not Gaussian, the knowledge of higher order cumulants is important to 
characterize the process. 

In the following we show that indeed the MFPT of Levy flights is far from being rep- 
resentative, and knowledge of higher order moments, or actually the whole distribution, is 
needed to give a proper description of the statistics of FPT. In order to achieve this task, 
we start by obtaining detailed knowledge about the spectrum of the fractional Laplacian. 
Using analytical and numerical results of the spectrum, we then show that higher order 
moments are important in comparison to MFPT, especially in the vicinity of the absorbing 
boundaries. We also obtain the whole distribution of FPT. Finally, we comment on how 
this approach can be generalized to other boundary conditions, such as reflecting or mixed 
ones. 

Fractional Laplacian in Bounded Domains. — In a recent paper [15] we developed 
an approach to study the spectrum of large powers of the Laplacian in bounded domains, 
continuing a previous effort [16] which focused on the ground state (i.e., the eigenfunction 
corresponding to the lowest eigenvalue). We showed that in the large TV-limit the eigen- 
functions of (— A) N with absorbing BC are simply proportional to the associated Legendre 
polynomials P2N+j( x ) f° r i ^ ^> an< ^ we showed how to obtain systematic corrections in 
powers of 1/7V. It turns out that this asymptotic expansion shows remarkable convergence, 
so that already for TV = 1 the asymptotic spectrum is very close to the exact one. In addi- 
tion, we showed that expressing (— A) N in the basis of the associated Legendre polynomials 
does not only diagonalize it for N — > oo, but is also a very useful basis for numerically 
evaluating its spectrum for any N. However, these results were limited to integer TV's and 
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are therefore inapplicable to the phenomena which interests us here, namely of anomalous 
diffusion. 

Here we show how to extend the method developed in [15] to non-integer N'a. The 
challenge is to find an appropriate basis that satisfies the BC (2), and coincides with the 
associated Legendre polynomials for integer TV's. The following normalized functions form 
the required basis 



_ r(4JV+l)y / (2j)! (27V+2j+±) (l-x 2 ) N (W+i) 

& ( X > = 4«r(2jv+i)Vr(4W+2 J+ i) °V W ' 

for j e N (including j = 0), and C^\x) arc the Gegenbauer polynomials [17]. In order to 
show this, we need to write the matrix elements of the operator (— A) N in this basis, namely 



£& d =J Jm (*) [(-A)" Si {*) dx = A£ m . (7) 

This can be achieved using -^^x x = r p^j_\) x x ~^ [7], and the series expansion of the 
Gegenbauer polynomials [17] 

G 2j £[(2j-2£)\ ' W 



1=0 



where (X) n is the Pochhammer symbol. After some algebra similar to [15], the following 
expression for the matrix elements are obtained 

S N / (2j)!(4jV+2m)!V(4jV+4m+l)(4W+4j+l)(-l) w +^ * (-4)' (2jV+2i)!r (2N+j+i+ 1) T (i+ 1) 

™ J V (2m)!(4JV+2j)! 2 ^ + ir(jV+i) £{, [(2i)!] 2 0-i)!r(iV+;+§) 



*+±,* + l 'V 32 



2' ■-)-■ . - 

Jtf,i + i y " " y zjv -i- l,N + i + | 

a, 6, c 



2iV + m + i,-m,iV+ 1 

3 
2 



where 3 F 2 ^ ^ ^ ;xj is the generalized hypergeometric function [17]. Interestingly, for 

integer N's this expression coincides with the one obtained in [15]. This implies immediately 
that all the results obtained for large N's apply also for the fractional case, and we can readily 
use the approach developed in [15] with the important difference that for non-integer TV's 
one should use as basis the functions (6) and not the associated Legendre polynomials. We 
therefore obtain for the cigcnfunctions 



[(2j)(2j-l)] 3/2 , 
Vj {x) = fj(x) + —— 2 fj-i{x) 



327V 2 
[(2 J + 2)(2j + l)] 3/2 
327V 2 



and for the eigenvalues 

= 



^w r(2iV+1) 



3 + 4j + 8,? 2 
167V 



75 + 344 j + 672j 2 + 832j 3 + 192j 4 
+ 15367V 2 + 



1 

TV 3 " 



(10) 



(11) 



In particular, the numerical scheme of [15] will allow us to calculate to arbitrary precision 
all the quantities of interest, and thus provides a reference to the analytical results we present 
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below. Since we are interested in anomalous diffusion and thus primarily on small values 
of N, we will focus from now on only on the case < N < 1. To demonstrate that the 
numerical scheme of [15] can be extended to noninteger as well as small N's, we present in 
Fig. 1 the first three eigenvalues Xo(N), Xi(N), X 2 (N) in the range < N < 1. As can be 
seen, the lowest eigenvalue Xq(N) becomes smaller than one for values of N < 0.32. 
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Fig. 1: The first three eigenvalues X (N), Xi(N) and X 2 (N) of (-A)^ in the range < N < 1. Note 
that all the eigenvalues are degenerate at N — where their value is 1, and that only Ao becomes 
smaller than one in a certain range of N (see inset). 



In the same line of [15] where a 1/iV expansion has been performed, we begin by studying 
the small N behaviour of the operator (— A) . To zeroth order in N, it is easy to show 
that Eq. (9) yields A^ j = 8 m j + O(N), which is consistent with the simple fact that 
(—A) = 1. Going beyond zeroth order is cumbersome, and is mainly due to the expansion 
of the Generalized Hypergeometric functions. However, we have previously shown that (see 
appendix B in Ref. [15]) 



2N + m+i,-m,N + l \ _ -A ?! B(2N+ m +t+\,\~i) B (N+i,i+i+\\ 

2N + 1, N + i + | ) ~ a(m - eV B{2N+m+\,\-m) B^N+l^j 4 ' 



i- j,-N,2N + + | \ _ 4- 2W -^- v ^(2i)!(j-z)! 



3^2 



j-i 

(iN+2j+2m+2i)\ JVJ 

/ j (2m+2i)! (j-m-i) ! (2 jV+ j +m+i) ! m\{N-m)\ 
m=0 



(13) 



where B(x,y) is the Beta function. These representations render the small N expansion 
straightforward yet tedious. Then, the expansion of A^ • up to first order in TV is given by 



with 



and 



imj - S mj + NB m j + O (N 2 ) , (14) 



j Dj m = j 

Bm,3 = { % /(4m+l)(4j+l) , . 7 (15) 

t |m(2m+l)-j(2j+l)| ' m r 3 



Dj = 4^ (4j) - 4^o (2j) - 7(4^T) - 2 7 

4.7+1 (- 4 )' r (j+'+|) ^ (-4)'r(<+3+|) \l-^+t+^)+^o{ 3 +l+\)} 
+ 2. L (2i)!(3-t)l (w)iy-^i ' [ ' 

i=0 £=0 V 2; 
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Table 1: The first eigenvalues of the matrix B given by Eq. (15) - denoted bj. Note that the lowest 
eigenvalue is negative (bo < 0), while all the rest are positive (bj>\ > 0). 



3 





1 


2 


3 


4 


5 


h 


-0.688 


2.727 


3.909 


4.646 


5.183 


5.606 



where ipo(x) is the digamma function and 7 is the Euler-Mascheroni constant [17]. Note 
that D is obtained by taking the limit j — > in Eq. (16), namely 

D a = 2(1 - 7 - In 2) ~ -0.541 . (17) 



Equations (15-17) allow determining the eigenvalues and the eigenfunctions of the frac- 
tional Laplacian for small N. Note that the zeroth order term in the expansion of (— A)^ is 
just the identity, implying that at the lowest order, the eigenvalues are completely degener- 
ate. This means that diagonalizing (— A) N to first order is equivalent to the diagonalization 

of the matrix B given by Eq. (15). Once the eigenvalues {bj} and eigenvectors {(3j} of the 
matrix B are computed numerically, the eigenvalues {A^} and eigenfunctions {vf(x)} of 
(— A) N for small N follows immediately and are just given by 



vf(x) 



1 + Nbi 



E/#°/»(*) : 

n=0 



(18) 
(19) 



where /3j is the n th component of f3j. Notice that at this order, the iV-dependence of v^(x) 
comes in only through the basis functions f n (x). 

In Table 1, we provide the first eigenvalues of B m j . Interestingly, only the smallest 
eigenvalue bo is negative, implying that A^ < 1, consistent with the numerical results 
presented in Fig. 1. The other eigenvalues bj>\ are all positive. This is in contradiction 
with the small N expression developed in [12] for the lowest eigenvalue, predicting Ao = 
1 + 2N(l - 7) > 1 (Eq. (36) there). The reason for this difference is twofold. First, in [12] 
the operator (— A) N has been written using a trigonometric basis yielding a different first- 
order matrix B m j. For example, in the present notation, the value of Do obtained by [12] 
is Dq RK = 2(1 — 7) ~ 0.846. Second, the matrix B m j has not been fully diagonalizcd 
in [12], but rather approximated by A = 1 + NDq HK , predicting a positive correction for 
Ao for small N's. This suggests that the trigonometric basis used in [12] is not adequate for 
a perturbative analysis. Effectively, if we use the same approximation Ao = 1 + NDq we 
would obtain a negative correction, implying that the basis {fj(x)} is more adapted to the 
operator (— A) N . Note that the fact that the lowest eigenvalue is smaller than one can have 
important consequences regarding the stability of the process described by the fractional 
Laplacian. To allow for a more advanced comparison, we also show in Fig. 2 the ground 
state Vq(x) obtained using this approximation, compared to the exact ground state v (x) 
computed numerically from the diagonalization of (9), for various small values of N. 

Complementary to the large- TV study of Ref. [15] and the small- study presented above, 
we now develop a useful approach for obtaining the ground state directly from the matrix 
Afj , without expanding around a particular value of N. Often in physical problems the 
ground state is the most interesting quantity, as will be shown to be the case for FPT too. 
For this purpose, we approximate the matrix A^ as given by Eq. (9) so that the lowest 
eigenvalue is given by diagonalizing its upper 2x2 sub-matrix. We then find for the ground 
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Fig. 2: A comparison of the exact ground state VQ Xact (x) with vq{x) for various low values of TV. 
The Insets show the difference VQ Xac *(x) - Vq(x). 



state 



where 



A^TV) 



r(7V)r(7V+ l)r(47V + 2) 
2 2N + 1 T(2N)T(2N + 2) 



ab , 



v 2 x2 (x) = [f (x)+af 1 (x)]/Vl + a\ 



b = 



(c + 1) - V(c-l) 2 +46 2 /26 , 
TV /47V + 5 



27V + 3 V 27V + 1 ' 

(47V + 5) (167V 4 + 407V 3 + 427V 2 + 187V + 3) 
(27V + l)(27V + 3)(27V + 5) 



(20) 
(21) 



(22) 



As can be seen in Figs. 3-4, these expressions compare well with the exact numerical diag- 
onalization of At! A . 
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Fig. 3: A comparison of the smallest eigenvalue \^ xact with A 2 ,* 2 resulting from the 2x2 approx- 
imation for various values of TV. Note that the relative error is always smaller than 2 • 10~ 4 (see 
Inset) and that A 2 ,* 2 becomes exact when TV — > oo [15]. 



First Passage Times. — Equipped with the new results concerning the spectrum of 
the fractional Laplacian, we can readdress the statistics of FPT of Levy flights. It is easy 
to solve formally Eq. (1) as 



C(x,t\x ) =J2 erXAN)tv 3 ( x "> v i 0"°) 

3=0 



(23) 
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where here too {vj(x)} are the eigenfunctions of (— A) . Using Eq. (3) we obtain 

P (<l*o) = £ ( / Vj ( x ) dx ) X i ( N ) vj (*o) e- X]{N)t , (24) 
i=o ^ J ~ 1 ' 

and Eq. (4) for the moments becomes 

00 ( J—i v i ( x ) dx) 
{t m ) (x ) =T(m + l)J2 K [A , {N)] m ' v, fro) . (25) 

We calculated numerically the moments using both Eq. (25) and a direct numerical 
diagonalization of (— A) N as given by Eq. (9). In Fig. 5, we present the cumulants of 
the distribution K m (xo) for N = 1/4,1/2,3/4. Recall that the cumulants simple 
combinations of (t m ) [17], and that for to = 1,2 they coincide with the mean and the 
variance respectively. We also present results obtained using only the first term in Eq. (25), 
namely the ground state, which we calculate here using the 2x2 approximation (20,21). As 
can be seen in Fig. 5 the agreement between the two approaches is very good. This means 
that knowledge of the ground state provides a reasonable description of the distribution. 
Note that the cumulants grow faster as the value of N decreases. 

However, these results show that Km™ increases with m for all values of xq, implying 
that the mean (t) = n\ is generically not representative, since the distribution is not peaked 
around it. In that case the variance, and other higher order cumulants, are necessary to 
characterize the full distribution. Put differently, statistical inference regarding first passage 
properties based only on the mean can lead to wrong conclusions since the errors bars, 
calculated using the higher order moments, are typically larger than the the mean value 
itself. Last, Fig. 6 shows the full distribution of FPT, p(t\xo), for N = 1/2, calculated using 
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Eq. (24) and a straightforward application of the numerical scheme presented above. For 
long times t 1 the tail of p(t\xo) is always exponential, and the lowest eigenvalue Ao(TV) 
controls its decay rate. However, for short times, many terms has to be retained in the sum 
(24) giving rise to a nontrivial form. 
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Fig. 5: A comparison of the first five cumulants as function of xo between the exact numerical result 
and the 2 x 2 approach for TV = 1/4, TV = 1/2 and TV = 3/4. The symbols ■, ♦, ▲, *, and • 

correspond to (t), 1 J A and k^ 5 respectively. 




Fig. 6: The distribution of FPT p(t\xa) as a function of xq and t for TV = 1/2. 



Summary and discussion. — In this letter we obtain both analytical and numerical 
information about the spectrum of the fractional Laplacian in bounded domains, extending 
previous results for integer powers of the Laplacian [15]. These results allow us to address 
the timely question of FPT statistics. We show that in general the mean first passage time 
does not provide a proper description of the FPT distribution since the variance (as well as 
higher order cumulants) become more and more important as the Levy stability index TV 
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becomes smaller (see Fig. 5). This means that the distribution is not peaked around the 
average, as often implied when only the average is discussed. Therefore, the MFPT is not 
representative of the distribution, and the whole distribution p(t\x a ) has to be retained. A 
direct consequence of these results is that the distribution of FPT is not Gaussian, contrarily 
to what is implicitly assumed when focussing on the average (MFPT) only. We show how to 
obtain the FPT distribution numerically for any N. In addition we show that a simple "2x2" 
approximation is able to reproduce accurately the cumulants. Our results (24) suggest that 
the tail of the distribution is exponential, while for short times it has a complex nontrivial 
form. 

An interesting point, is that the lowest eigenvalue A becomes smaller than one (see the 
inset in Fig. 1) in a whole range of small N's. Even though this result does not affect the 
FPT statistics studied here, this might have important consequences in systems described 
by the Fractional Laplacian, whose stability is determined by large powers of the eigenvalue 
of the ground state (an example of such a system can be found in [18]). Our result suggests 
that in such cases tuning the value of N (by changing the fractality of the medium for 
example) can control the stability of the system. 

Note that we were interested above only in the even part of the spectrum of (— A)^ 
(that is the even eigenfunctions). The reason is that for the statistics of FPT with two 
absorbing boundaries we need to know only the even spectrum since the odd eigenfunction 
do not contribute to the expression for p(t\xo) as given by Eq. (24) due to the integral 
there. However, information about the odd spectrum can be easily obtained by replacing 



Last, we comment on how to treat other boundary conditions, such as reflecting or mixed 



ones [11]. The key observation is that one has C> 2 (x) — Pj ' (x) where P^' v ' (x) is 



the j th Jacobi polynomial [17]. Using orthogonality properties and BC of this polynomials 
it can be shown that the relevant basis for mixed BC (i.e., one reflecting and one absorbing 
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